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ABSTRACT 

We develop a new Bayesian method for estimating white noise levels in CMB sky maps, and apply 
this algorithm to the 5-year WMAP data. We assume that the amplitude of the noise RMS is scaled by 
a constant value, a, relative to a pre-specified noise level. We then derive the corresponding conditional 
density, P{a \ s, Cg, d), which is subsequently integrated into a general CMB Gibbs sampler. We first 
verify our code by analyzing simulated data sets, and then apply the framework to the WMAP data. 
For the foreground-reduced 5-year WMAP sky maps and the nominal noise levels initially provided 
in the 5-year data release, we find that the posterior means typically range between a = 1.005 ±0.001 
and a = l.OlOiO.OOl depending on differencing assembly, indicating that the noise level of these maps 
are biased low by 0.5-1.0%. The same problem is not observed for the uncorrected WMAP sky maps. 
After the preprint version of this letter appeared on astro-ph., the WMAP team has corrected the 
values presented on their web page, noting that the initially provided values were in fact estimates from 
the 3-year data release, not from the 5-year estimates. However, internally in their 5-year analysis the 
correct noise values were used, and no cosmological results are therefore compromised by this error. 
Thus, our method has already been demonstrated in practice to be both useful and accurate. 
Subject headings: cosmic microwave background — cosmology: observations — methods: numerical 



1. INTRODUCTION 

The cosmic microwave background (CMB) is probably 
the most valuable source of observational data in modern 
cosmology. Several experiments have been carried out to 
map its anisotropics, most notably the Wilkinson Mi- 
crowave Aniso tropy Map (WMAP) (jBennett et al.ll2003l : 
iHinshaw et~ani2007ft . The WMAP experiment has pro- 
vided unique new insights in the workings of the universe, 
from large to small scales, and we now believe that we 
understand the main physical process from after inflation 
and up until today. 

The theory of inflation was initially propo sed as a so- 
lution to the horizon and flatness problem (jGuth et al I 
Il981f ). Additionally, it established a highly suc- 
cessful theory for the formation of primordial den- 
sity perturbations, thus providing the required seeds 
for the large-scale structures (LSS), later giving rise 
to the temperature anisotropies in the cosmic mi- 
crowave backgrou n d radiation that we observe today 
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From this theory, we are able to predict what the sta- 
tistical properties of the CMB map should be, given a set 
of cosmological parameters. The goal of the cosmological 
data analyst is then to determine how well this universe 
model fits real-life data, which is contaminated by fore- 
grounds, systematics and various uncertainties. The end 
result may be summarized in terms of a joint posterior 
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including all unknown quantities, from which the desired 
cosmological parameters may be obtained by marginal- 
izing over any relevant nuisance parameters. 

A long-time discussion within the field of cosmolog- 
ical data analysis has revolved around the process of 
optimally extracting the underlying cosmological sig- 
nal from the data. One general and potentially op- 
timal framework for doing this, based on a statis- 
tical algo rithm called Gibb s sampling, was first de- 
scribed byEewell et all (|2004D . I Wandelt et al.l(|2004f ) and 
Eriks en et al.l ( 2004|). and extend e d and applied to the 
WMAP d ata by lO 'Dwver ct al.' ("2004'); 'Eriksen et aL| 
(l2bG7a,b,) : ILarson et al. ( _2007) ; Eriksen et al. (2008a, 1^; 
iGroeneboom fc EriksenI ( 20091 ). This algorithm provides 
the user with samples drawn from a joint CMB poste- 
rior, which can also include a large number of external 
nuisance parameters. One well-known and important ex- 
ample of this is joint component separation and CMB 
power spectrum estimation. 

We have developed an independent implementation of 
the CMB Gibbs sampler which we call "Slave" , corre- 
sponding to a light weight C -I- -I- v ersion of "Comman- 
der", described bv lEriksen et al.l ()2004f ). While Slave 
may lack advanced features such as foreground estima- 
tion and multi-band analysis, it does include all the basic 
features required for elementary CMB power spectrum 
analysis (e.g., support for cut-sky analysis, anisotropic 
noise distribution) , and it also takes advantage of the ob- 
ject oriented design features available in C-I--I-. This is 
particularly useful when adding additional features into 
the Gibbs sampler, as exemplified in the present paper. 

In this paper, we apply this new implementation to a 
new practical problem, and consider estimation of the 
white noise level of CMB sky maps directly from the 
data. Specifically, we develop the necessary machinery 
for including this operation into the CMB Gibbs sampler, 
and apply this tool to the 5-year WMAP data. 
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Fig. 1. — Posterior distribution P{a\d) for the foreground- 
corrected W3 channel. This channel exhibits the strongest de- 
tection of a ^ 1 of any DA. 

2. METHODS 

An observed CMB map may be modeled as: 

d = As + n (1) 

where A denotes the instrumental beam, s is the desired 
CMB signal, and n is instrumental noise. The noise is 
in this paper assumed to be uncorrclatcd, and the corre- 
sponding noise covariance matrix in pixel space is thus 
Nij = afSij, where ai is the noise standard deviation for 
the ith pixel. Further, we assume that the CMB fluctua- 
tions are Gaussian and isotropic, so the signal covariance 
matrix simplifies to Cem.i'm' = CeSu>6mm'- 

In order to estimate the power spectrum Cg and the 
signal s given the data, we need to sample from the 
joint distribution P{Ci, s\d). The algorithm for sampling 
from this dis tribution by Gibbs sampling is describe d 
exstensivly bv lJewell et"al] ll2004l).IWandelt et all (|2004f ) , 
citctchu:2005 and lEriksen et all (f2004f ). A self-contained 
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the Gibbs sampling algo- 
(|2009( 1. together with 



rithm is presented in iGroeneboor 
a presentation of the "Slave" framework. 

The standard Gibbs sampler draws samples from 
the joint distribution, P{s, Ct\d), by alternately sam- 
pling from the conditional distributions P{s\C'i,d) and 
P{Ci\s). If one wishes to introduce further parameters 
into the data model, all that is required for joint estima- 
tion with the existing parameters is a sampling algorithm 
for the corresponding conditional distribution. 

Traditionall y, the noise properti es used in the Gibbs 
sampler (e.g., lEriksen et al.l [200^ have been assumed 
known to infinite precision. In this paper, however, we 
relax this assumption, and introduce a new free param- 
eter, a, that scales the fiducial noise covariance matrix, 
A'''''^,, such that N = aN^'^. Thus, if there is no devia- 
tion between the assumed and real noise levels, then a 
should equal 1. 

The full joint posterior, P(s,Ci,a\d), now includes 
the amplitude a. We can rewrite this as follows: 

P(s, Ci, a\d) = Pid\ s, a) ■ P{s, d) ■ P{a) (2) 

where the first term is the likelihood. 



P{d\s,a) 
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(3) 



the second term is a CMB prior, and the third term is 
a prior on a. Note that the latter two arc independent, 
given that these describe two a-priori independent ob- 
jects. In this paper, we adopt a Gaussian prior centered 



on unity on a, P{a) ^ N{1,(t'^). Typically, we choose 
a very loose prior, such that the posterior is completely 
data-driven. 

The conditional distribution for a can now be ex- 
pressed as 



P{a I s, Ci, d) (X 



i/2 



P{a) 



(4) 



where n = Npi^ and (3 = {d — s)N ^{d — s) is the x^. 
(Note that the is already calculated within the Gibbs 
sampler, as it is used to validate that the input noise 
maps and beams are within a correct range for each 
Gibbs iteration. Sampling from this distribution within 
the Gibbs sampler represent therefore a completely neg- 
ligible extra computational cost.) For the Gaussian prior 
with unity mean and standard deviation ctq, , we find that 



P{a I s, C^, d) cx 



yn/2 



(5) 



For large degrees of freedom, n, the inverse gamma 
function converges to a Gaussian distribution with mean 
fi ~ b/{k + 1), where we have defined k — n-p\y^/2 — 1, 
and variance a"^ — b'^/{{k — \){k — l)(fc — 2)). A good 
approximation is therefore letting ai+i be drawn from a 
product of two Gaussian distributions, which itself is a 
Gaussian, with mean and standard deviation 
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(6) 



(7) 



This sampling step has been implemented in "Slave" , 
and we have successfully tested it on simulated maps. 
With A"side = 512 and ^max = 1300 and full sky coverage, 
we find a = l.OOOzbO.OOl. Note that with such high reso- 
lution, the standard deviation on a is extremely low, and 
any deviation from the exact a = 1.0 will be detected. 

3. DATA 

In this paper we an alyze the 5-year WMAP 
data ()Hinshaw et"aI1 l2009f ). which is available from 
LAMBDA". We consider both the raw sky maps for 
each differencing assembly (D A), and the corre spond- 
ing foreground-reduced maps (IGold et all [2001 . The 
nominal noise amplitudes are taken from the original 5- 
year data release as presented on LAMBDA. However, 
we note that after the initial publication of this paper, 
these have been corrected, as there was a discrepancy be- 
tween the values presented on LAMBDA and those pub- 
lished in the WMA P Five Year Explanatory Supplement 
(jLimon et al.ll200"9[ l. The updated values agree very well 
with our results. 

The foreground-reduced maps were produced from the 
raw maps by fitting and subtracting three fixed tem- 
plates to each case. One of these templates was the 
(K-Ka) difference map, which mainly traces synchrotron, 
free-free and, possibly, spinning d ust, and the other two 
were the Ha template of lFinkbeincr ( 2003) and the FDS8 
thermal dust template of iFinkbeiner et all (|1999( ). Here 

http://lambda.gsfc.nasa.gov 
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TABLE 1 
Noise estimation results 





Oi — 


1 (in %) 




(mK) 




Band 


Direct 


Bayesian 


Nominal 


Estimated 


Updated 




Raw maps 


Kl 


1.90 


0.50 


1.436 


1.439 




Kal 


0.10 


0.10 


1.470 


1.471 




Ql 


0.23 


0.15 


2.254 


2.256 




Q2 


0.16 


-0.10 


2.141 


2.140 




VI 


0.08 


-0.08 


3.314 


3.313 




V2 


-0.06 


-0.03 


2.953 


2.953 




Wl 


0.52 


-0.37 


5.899 


5.889 




W2 


0.23 


-0.43 


6.565 


6.550 




W3 


0.33 


-0.30 


6.926 


6.916 




W4 


-0.45 


-0.85 


6.761 


6.732 




Foreground-reduced sky maps 


Ql 


0.83 


0.51 


2.2449 


2.2542 


2.254 


Q2 


0.87 


1.02 


2.1347 


2.1455 


2.141 


VI 


0.71 


0.50 


3.3040 


3.3123 


3.314 


V2 


0.63 


0.62 


2.9458 


2.9549 


2.953 


Wl 


1.06 


0.53 


5.8833 


5.8988 


5.899 


W2 


1.27 


0.71 


6.5324 


6.5555 


6.565 


W3 


1.54 


1.12 


6.8849 


6.9233 


6.926 


W4 


0.04 


-0.20 


6.7441 


6.7373 


6.761 



Note. — Summary of 5-year WMAP noise amplitudes. Sec- 
ond column: Noise amplitudes estimated directly from high-£ 
power spectra, quoted in terms of o — 1 in percent. Third 
column: Noise amplitudes estimated by Gibbs sampling. The 
uncertainty on each of these numbers is 0.1%. Fourth column: 
Noise RMS per observation as quoted by the WMAP team 
on LAMBDA. Fifth column: Noise RMS per observation es- 
timated by Gibbs sampling, sixth column: Noise RMS per 
observation as quoted by the WMAP team in the Five Year 
Explanatory Supplement, page 65. 

it is worth noting that the (K-Ka) diffe rence map was 
smo othed to 1° FWHM before fitting (jHinshaw et al.l 
120071) ; ahhough this difference map is intrinsically noisy, 
it does not have power beyond I ^ 400, and, in effect, 
the raw and foreground-corrected maps are identical at 
high £'s. This will be explicitly demonstrated later. 

In the low signal-to-noise regime, the estimated CMB 
signal will fluctuate greatly. In order to dampen the un- 
ruly behavior in this regime, we have chosen to bin the 
power spectrum for high £s. The C^s are then generated 
from the binned signal power spectrum. 

Unless explicitly noted, we apply the KQ85 WMAP 
sky cut (jGold et al.1 l2009f ) . which removes 18% of the 
sky, including point source cuts. We also take into ac- 
count the circular-symmetric beam profiles for each DA. 
Finally, the main analysis is carried out at a HEALPix^ 
resolution of iVgido = 512 and included harmonic space 
multipolcs between ^ = 2 and 1300. Including such high 
multipoles in the multipolc expansion is acceptable in 
this case for two reasons; first, the WMAP beams fall 
off quickly in this regime, and the data becomes strongly 
noise dominated. Second, we bin the angular power spec- 
trum heavily at high 

4. RESULTS 

The main results from our analysis are given in the 
third column of Tabic [TJ where the posterior mean of 
a is given for each DA, both for raw and foreground- 
corrected maps. The numbers are quoted in terms of 

^ http://hcalpix.jpl.nasa.gov 



a — 1 in percent, and the corresponding posterior RMS 
is 0.1%. 

For most of the raw DA sky maps, we find a = 1 
to within a few sigma. The largest outlier is W4, with 
a negative amplitude of -0.85%. This DA is known to 
have the strongest correlated noise of any WMAP DA. 
In general, we find that the noise levels of the raw sky 
maps are in good agreement with the levels quoted by 
the WMAP team. 

However, the situation is different for the foreground- 
corrected maps relative to the original noise amplitudes 
presented on LAMBDA. Specifically, in general the am- 
pfitudes of these maps are shifted high by 0.5 — 1.0%. 
The most extreme shift is seen for the W3 DA, with a 
posterior distribution centered on a = 1.011 ± 0.001. 

With such significant discrepancies between the pre- 
dicted and observed noise levels for the foreground- 
reduced sky maps, it should be possible to observe an 
absolute shift in the high-^ power spectra between simu- 
lated and the real sky maps. We therefore implemented 
a simple and approximate method to estimate a directly, 
without going through the Gibbs sampler: First, we cal- 
culate the angular power spectrum for each DA to high 
£'s, where the beam has killed all signal, and only noise 
is left. Any sky cut is ignored in this step. We then 
compute the average spectrum amplitude at ^ > 1300, 
and define this as our noise amplitude estimate. Next, 
we simulate an ensemble of noise realizations from the 
RMS maps of each DA, and repeat the above calculation 
for each of these, and compute the corresponding aver- 
age. The ratio between observed noise spectrum ampli- 
tude and the simulated average is then an estimate of 
a, and can be compared to the values obtained from the 
Bayesian analysis. Again, note that this is only a rough 
cross-check on the our results, and should not be consid- 
ered a proper stand-alone result because the sky cut is 
completely ignored. 

The results from this exercise are tabulated in the sec- 
ond column of Table[TJ In general, we see that these gen- 
erally follow the same trends as those from the Bayesian 
analysis, although with systematically slightly higher 
amplitudes. This is likely due to neglecting of the galaxy 
cut, as foregrounds may obviously add to the total spec- 
trum level. 

In figure [2] we plot the power spectrum of the raw 5- 
ycar WMAP W3 sky map in Ci units, together with the 
average W3 RMS cto value, as given by the WMAP team. 
Wc also plot the estimated new W3 RMS, as estimated 
by Slave. Here it is clear that the Slave estimate fits the 
power spectrum better. 

Finally, in Figure [3] wc plot the same information, but 
to lower £'s and on a logarithmic scale. Notice how 
the high-£ W3 foreground reduced maps power spectrum 
equals the spectrum from the raw map. This implies 
that both the foreground-reduced and raw data have the 
same high-€ noise levels, and do not differ. We therefore 
conclude that the shift in the noise amplitude is not due 
to any shift in the actual data, but rather a difference in 
the WMAP5 noise models. 

In the fourth column of Table [H we tabulate for easy 
reference the noise RMS values for a single observations 
listed on LAMBDA. Here it is seen that the predicted 
noise levels for the foreground corrected maps is indeed 
lower than those for the raw maps. However, this does 
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WMAP W3 raw/reduced 

WMAP W3 RMS 

Slave estimated W3 RMS 




Fig. 2. — The raw/foreground reduced power spectrum (red) 
plotted with the WMAP noise (green) and the Slave-estimated 
noise (blue). Note how Slave estimates a different noise level (blue) 
than the input noise level (green) . 



W3 raw data 

W3 noise power spectrum 
W3 foreground-reduced data 
Raw W3 noise power spectrum 




Multipole I 

Fig. 3. — A noise power spectrum realization for W3 RMS (black) 
and raw W3 RMS (red). The data power spectrum can be seen 
entering at high £ (blue and green) , where the raw and foreground- 
corrected data converge. 

not seem to be the case, neither from looking at the raw 
power spectra, nor from the results from our Bayesian 
analysis. 

Again, following the first publication of this letter, it 
was found that the foreground-correct noise RMS val- 
ues on the LAMBDA website were not properly updated 
with the 5-year WMAP results. Therefore, the noise lev- 
els used in this analysis are known to be incorrect, and 
the reason is fully understood. Thus, this demonstrates 
the new method works well as the estimated values are 
fully consistent with the corrected 5-year values, shown 
in the sixth column of Table [1] column 6. 

5. CONCLUSION 

We have introduced a new method for estimating white 
noise levels in CMB sky maps, using a Bayesian frame- 
work. We have then applied this method to the 5-year 



WMAP data, and re-estimated the noise levels of both 
the raw and foreground-reduced sky maps. In doing 
so, we found that the predicted noise levels for the raw 
maps are in acceptable agreement with the predictions, 
while the noise levels in the foreground-reduced maps are 
0.5 — 1.0% higher than the estimate initially provided by 
the WMAP team on LAMBDA. 

The explanation for this effect has after the publica- 
tion of this paper been found by the WMAP team simply 
to be an error in the results provided on LAMBDA: The 
quoted values were derived from the 3-year analysis in- 
stead of the 5-year analysis. However, the correct values 
were used in their cosmological analysis for the 5-year 
data, and no results are therefore compromised by this 
error. 

Thus, the method has already been demonstrated to 
be both accurate and useful on a practical example. Fur- 
ther, it carries virtually no extra computational cost 
within a Gibbs sampler, since all required quantities are 
already computed within this algorithm. We therefore 
recommend this feature to be used as a standard part 
of the Gibbs sampling machinery, since it both provides 
additional robustness against noise mis-estimation, and 
also propagates the uncertainty in these estimates into 
the final CMB power spectrum. 

Finally, we note that these noise estimates are very ro- 
bust with respect to systematic issues such as foreground 
or CMB signal estimation. The reason is that the sky 
maps are well oversampled; with ^max — 1300 there are 
~ 1.7 million modes in the harmonic expansion, while at 
A'sidc = 512 there arc ~ 3 million pixels. Therefore, there 
are essentially 1.3 million modes available to estimate one 
single number, a. If even greater precision is desired, one 
could simply consider increasing the pixel resolution one 
step further, which would leave the sky signal unchanged, 
because it is bandwidth limited, whereas the total num- 
ber of modes in the map quadruples, thus decreasing the 
noise uncertainty. 
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